a=0; # da ne dojame kot funkcijo

# ----------------------------------------------------
# izracunaj utez
# ----------------------------------------------------
function [w]=utez(x)
  w=-6*x^5 + 15*x^4 - 10*x^3 + 0*x^2 + 0*x^1 + 1;
  #w=1-x;
endfunction

# ----------------------------------------------------
# ustvari matriko visin za GLAVNE TOCKE, ter matriki za X in Y 
# koordinate za vsako GLAVNO tocko,
# velikost je (stGlavnihKvY+3)*(stGlavnihKvX+3)
# ----------------------------------------------------
function [matrikaVisin,matrikaPozicijX,matrikaPozicijY]=napolniMatrikoGlavnihTock(stGlavnihKvadratovX,minX,maxX)
  velikostGlavnegaKvX=(maxX-minX)/stGlavnihKvadratovX;
  stTockX=stGlavnihKvadratovX+1;
  matrika=zeros(stTockX,stTockX);
  offsetX=-(maxX-minX)/2;
  currentX=offsetX;
  currentY=offsetX;
  for i=[1:stTockX]
    for j=[1:stTockX]
      # tukaj lahko vnesemo poljubno funkcijo ce zelimo (zakomentiraj rand(..) spodaj)
      #matrikaVisin(i,j)=enacbaPovrsine(currentX,currentY);
      matrikaPozicijX(i,j)=currentX;
      matrikaPozicijY(i,j)=currentY;
      currentX+=velikostGlavnegaKvX;
    endfor
    currentX=offsetX;
    currentY+=velikostGlavnegaKvX;
  endfor  
  matrikaVisin=rand(stTockX,stTockX);
endfunction

# ----------------------------------------------------
# narisi mrezo
# ----------------------------------------------------
function narisi(mreza,stKvX,minX,maxX)
  tx = linspace (minX, maxX,stKvX+1);
  ty = linspace (minX, maxX,stKvX+1);
  [xx, yy] = meshgrid (tx, ty);
  mesh (tx, ty, mreza);
endfunction

# ----------------------------------------------------
# izracunaj odvod za dano tocko (za tocke glavnih kvadrator)
# ----------------------------------------------------
function [u_x,u_y,u_xx,u_yy,u_xy]=izracunajOdvodeZaDanoTocko(mrezaGlavnihKvadrator,indexX,indexY,velGlKvX)
  sz=size(mrezaGlavnihKvadrator);  
  m=mrezaGlavnihKvadrator;
  h=velGlKvX;
  yIndexFac=1;
  
  if((indexX+1)>sz(1,2))
    indexX=sz(1,2)-1;
  endif
  
  if((indexX-1)<1)
    indexX=2;
  endif
  
  if((indexY+1)>sz(1,1))
    indexY=sz(1,1)-1;
  endif
  
  if((indexY-1)<1)
    indexY=2;
  endif
  
  u_x=m(indexY,indexX+1)-m(indexY,indexX-1);
  u_x=u_x/(2*h);
  
  u_y=m(indexY+1*yIndexFac,indexX)-m(indexY-1*yIndexFac,indexX);
  u_y=u_y/(2*h);
  
  u_xx=m(indexY,indexX+1)+m(indexY,indexX-1)-2*m(indexY,indexX);
  u_xx=u_xx/(h*h);
  
  u_yy=m(indexY+1*yIndexFac,indexX)+m(indexY-1*yIndexFac,indexX)-2*m(indexY,indexX);
  u_yy=u_yy/(h*h);
  
  u_xy=m(indexY+1*yIndexFac,indexX+1)-m(indexY+1*yIndexFac,indexX-1)-m(indexY-1*yIndexFac,indexX+1)+m(indexY-1*yIndexFac,indexX-1);
  u_xy=u_xy/(h*h);
  
endfunction

# ----------------------------------------------------
# funkcijske vrednosti za f00,f10,f01,f11
# normX in normY sta [0,1] x in y pa dejanske koordinate na robu
# ----------------------------------------------------
function [v]=f00(x,y,z,xv,yv,matrikaVisin,xi,yi,velKv)

  [dx,dy,ddx,ddy,dxy]=izracunajOdvodeZaDanoTocko(matrikaVisin,xi,yi,velKv);

  a=ddx/2;
  b=ddy/2;
  c=dxy;
  d=dx;
  e=dy;
  v=a*x^2+b*y^2+c*x*y+d*x+e*y+z;

endfunction
function [v]=f10(x,y,z,xv,yv,matrikaVisin,xi,yi,velKv)
  x=1-x;
  
  [dx,dy,ddx,ddy,dxy]=izracunajOdvodeZaDanoTocko(matrikaVisin,xi,yi,velKv);  
  
  a=ddx/2;
  b=ddy/2;
  c=-dxy;
  d=-dx;
  e=dy;
  v=a*x^2+b*y^2+c*x*y+d*x+e*y+z;
  
endfunction
function [v]=f01(x,y,z,xv,yv,matrikaVisin,xi,yi,velKv)
  y=1-y;
  
  [dx,dy,ddx,ddy,dxy]=izracunajOdvodeZaDanoTocko(matrikaVisin,xi,yi,velKv);

  a=ddx/2;
  b=ddy/2;
  c=-dxy;
  d=dx;
  e=-dy;
  v=a*x^2+b*y^2+c*x*y+d*x+e*y+z;
  
endfunction
function [v]=f11(x,y,z,xv,yv,matrikaVisin,xi,yi,velKv)
  x=1-x;
  y=1-y;
  
  [dx,dy,ddx,ddy,dxy]=izracunajOdvodeZaDanoTocko(matrikaVisin,xi,yi,velKv);

  a=ddx/2;
  b=ddy/2;
  c=dxy;
  d=-dx;
  e=-dy;
  v=a*x^2+b*y^2+c*x*y+d*x+e*y+z;
  
endfunction
# ----------------------------------------------------
# interpoliraj
# ----------------------------------------------------
function [matrikaVsehTock]=interpoliraj(matrikaVisin,matrikaXCoord,matrikaYCoord,stGlKv,stPodKv,velikostGlKv)

  stVsehKv=stGlKv*stPodKv;
  matrikaVsehTock=zeros(stVsehKv+1,stVsehKv+1);
  size(matrikaVsehTock)
  # za vsak glavni kvadrat Y,X
  for gi=[1:stGlKv]
    for gj=[1:stGlKv]
    
      # iXY
      # | 00   01 |
      # | 10   11 |
      
      # indexi za v matriko glavnih tock (NE UPORABI ZA V PODTOCKE ZARADI ZAMIKA!)
      # giXX[x,y]
      gi00=[gj,gi];
      gi01=[gj,gi+1];
      gi10=[gj+1,gi];
      gi11=[gj+1,gi+1];
      
      x00=matrikaXCoord(gi00(1,2),gi00(1,1));  
      y00=matrikaYCoord(gi00(1,2),gi00(1,1));  
      x10=matrikaXCoord(gi10(1,2),gi10(1,1));  
      y10=matrikaYCoord(gi10(1,2),gi10(1,1));  
      x01=matrikaXCoord(gi01(1,2),gi01(1,1));  
      y01=matrikaYCoord(gi01(1,2),gi01(1,1));  
      x11=matrikaXCoord(gi11(1,2),gi11(1,1));  
      y11=matrikaYCoord(gi11(1,2),gi11(1,1));  
      
      v00=matrikaVisin(gi00(1,2),gi00(1,1));
      v01=matrikaVisin(gi01(1,2),gi01(1,1));
      v10=matrikaVisin(gi10(1,2),gi10(1,1));
      v11=matrikaVisin(gi11(1,2),gi11(1,1));
      
      
      
      stepSize=1/stPodKv;
      
      # normalized current values [0,1] 
      yParam=0;
      xParam=0;
      

      # za vsak podkvadrat Y,X
      for pi=[1:stPodKv+1]
        for pj=[1:stPodKv+1]
        

        
          # indexi za v matriko podtock (podkvadratov)
          # piXX[x,y]
          jj=gj-1;
          ii=gi-1;
          pi00=[jj*stPodKv+pj+0,ii*stPodKv+pi+0];
          pi01=[jj*stPodKv+pj+0,ii*stPodKv+pi+1];
          pi10=[jj*stPodKv+pj+1,ii*stPodKv+pi+0];
          pi11=[jj*stPodKv+pj+1,ii*stPodKv+pi+1];
                  
          w00x=utez(xParam);
          w00y=utez(yParam);
          
          w10x=utez(1-xParam);
          w10y=utez(yParam);

          w01x=utez(xParam);
          w01y=utez(1-yParam);
          
          w11x=utez(1-xParam);
          w11y=utez(1-yParam);
        
          
          interpX=(1-xParam)*x00+xParam*x10;
          interpY=(1-yParam)*y00+yParam*y01;
          

          u=f00(xParam,yParam,v00,x00,y00,matrikaVisin,gi00(1,1),gi00(1,2),velikostGlKv)*w00x*w00y+ f01(xParam,yParam,v01,x01,y01,matrikaVisin,gi01(1,1),gi01(1,2),velikostGlKv)*w01x*w01y+ f10(xParam,yParam,v10,x10,y10,matrikaVisin,gi10(1,1),gi10(1,2),velikostGlKv)*w10x*w10y+ f11(xParam,yParam,v11,x11,y11,matrikaVisin,gi11(1,1),gi11(1,2),velikostGlKv)*w11x*w11y;

          matrikaVsehTock(pi00(1,2),pi00(1,1))=u;

          xParam+=stepSize;
        endfor
        xParam=0;
        yParam+=stepSize;
      endfor
      
      
      
    endfor
  endfor
  
endfunction
# ----------------------------------------------------
# main
# ----------------------------------------------------

# MIN/MAX X IN Y MORATA BITI ENAKA, ST. X IN Y (POD)KV TUDI ENAKA

# ----------------------
#uporabnikovi parametri
# ----------------------
minimumX=-3; # zacetek X & Y
steviloKvadratovVSmeriX=4;
steviloPodkvadratorEnegaKvadrataVSmeriX=10; # 1 je isto kot "brez podkvadratov"


maximumX=-minimumX; # konec X
minimumY=minimumX; # zacetek Y
maximumY=maximumX; # konec Y
stVsehKvadratovVSmeriX=steviloKvadratovVSmeriX*steviloPodkvadratorEnegaKvadrataVSmeriX;

# ----------------------
# napolni matriko s vrednostmi visin, x coordinat in y coordinat
# (3 matrike, v vsaki ena komponenta (x,y,z)  za vsako tocko
# ----------------------

[matrikaVisinGlavnihTock,matrikaXCoord,matrikaYCoord]=napolniMatrikoGlavnihTock(steviloKvadratovVSmeriX,minimumX,maximumX);

# ----------------------
# izracunaj interpolacijo
# ----------------------

velikostGlavnegaKv=(maximumX-minimumX)/steviloKvadratovVSmeriX;

matrikaVsehTock=interpoliraj(matrikaVisinGlavnihTock,matrikaXCoord,matrikaYCoord,steviloKvadratovVSmeriX,steviloPodkvadratorEnegaKvadrataVSmeriX,velikostGlavnegaKv);


# ----------------------
# narisi
# ----------------------
#narisi(matrikaVisinGlavnihTock,steviloKvadratovVSmeriX,minimumX,maximumX);
narisi(matrikaVsehTock,stVsehKvadratovVSmeriX,minimumX,maximumX);
